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CRACKED FINITE ELEMENTS PROPOSED FOR NASTRAN 

By J. A. Aberson, The Lockheed-Georgla Company 
and J. M. Anderson, Georgia Institute of Technology 


SUMMARY 


The recent introduction of special crack-tip singularity elements, 
usually referred to as cracked elements, has brought the power and flexibility 
of the finite -element method to bear much more effectively on fracture mechan- 
ics problems. This paper recalls the development of two cracked elements and 
presents the results of some applications proving their accuracy and economy. 
Judging from the available literature on numerical methods in fracture mechan- 
ics, it seems clear that the elements described have been used more extensive- 
ly than any others in practical fracture mechanics applications. 


INTRODUCTION 


The study of crack growth behavior by classical continuum linear fracture 
mechanics has been limited primarily to simple structural configurations and 
loadings. Because of the ease with which the finite-element method handles 
discontinuous loads and boundary conditions, attempts have been made to use 
this method to study fracture in complex structures. The capability of con- 
ventional elements and conventional modeling procedures to predict crack-tip 
parameters accurately has been found to be limited and uneconomical. The most 
significant and accurate results obtained through use of standard methods are 
those which follow the work reported by Chan et. al. (ref. 1) and Kobayashi 
et. al. (ref. 2). The method followed by Kobayashi, for example, is to use the 
crack -opening displacements to solve a sequence of two simultaneous equations 
relating displacements to the two stress-intensity factors, Kj and Kn,for the 
opening and sliding fracture modes, respectively. These results are then in- 
terpreted or extrapolated to predict values at the crack tip. The difficulty 
in applying this approach lies in trying to accurately depict the extreme 
stress gradient existing in the near vicinity of a crack tip. The detail re- 
quired for a model to reasonably characterize this gradient makes the proce- 
dure expensive and cumbersome. For example, Oglesby and Lomacky (ref. 3) in- 
dicate that the maximum permissible element size necessary to Insure acceptable 
accuracy in the computed stress-intensity factors is of the order of l/300th 
of the crack half-length. To achieve detail of this order, substructure 
analyses have usually been used. That is, a coarse model is first analyzed to 
obtain boundary conditions which are then imposed on a more refined local 
model of the crack region. In some instances results from the second model 
have been used to analyze a third and even more finely modeled localized region 
of the crack tip. Obviously, considerable modeling and computer efforts are 
necessary to carry out such analyses. 
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Another approach, which is similar to that us^d by Chan and Kobayashi, is 
a combination of finite-element and boundary-collocation analyses. In this 
method, displacements and stresses from a conventional finite element model - 
which may or may not contain some representation of a crack - are used as 
boundary points in a boundary collocation solution for the crack. This is the 
procedure followed by Freese and Kaldjian (ref. 4), for example. The disadvan- 
tage of this combined approach is that it is not easily applicable since con- 
siderable experience in fracture analysis and complex variables is required to 
obtain consistently accurate results. 

In an attempt to circumvent the economic problem of the conventional 
approach and the applicability problem of the combined approach, research and 
development efforts have been turned toward formulating elements which are 
capable of characterizing the crack-tip stress singularity internally. These 
singularity elements - cracked finite elements - provide a new means for com- 
puting stress -intensity factors and thereby predicting crack growth. Two such 
elements have been developed and implemented at the Lockheed-Georgia Company. 
These elements have received extensive usage in project and contract work and 
have provided considerable substantiation of the accuracy of the approach used 
in their formulation. A description of the formulation and implementation of 
these elements follows. 


TECHNICAL DEVELOPMENT 


The Williams* series of stress functions (refs. 5 and 6)* is the basis 
for all boundary collocation and cracked finite-element schemes to estimate 
stress-intensity factors. This series gives the following familiar expressions 
for the in-plane stresses for the plane crack problem illustrated in figure 1. 
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* An error in reference 5 was subsequently corrected by Williams in 
reference 6. 


5>2 



and 


00 


T te ^,e) 


— -1 

n 2 1 f 


/n 


/ 4 r \ s n ' L < n+2 ) s in + 1; 9 

n=l 


f(n)(n-2) sin - lj0 


r fn 

+ a R g(n)(n+2) cos +1)0 - (n-2) cos 


(! • 1 >]j 


in which 


f(n) = 


t * 1 

(-D n 


and 


( 2 ) 


g(n) = 


f - (-l) n 

a + 1 

2 


The independent constants associated with the symmetric (even in 0) and anti- 
symmetric (odd in 0) parts of the Williams' series have been denoted in equa- 
tion (1) by s r and a R , respectively. Even though an expression for the second 

antisymmetric term has been formally written, it should be noted before pro- 
ceeding that it never contributes to any of the stresses. Thus in the usual 
finite-element description, a ^ is not a legitimate generalized coordinate. 

- 1/2 

The leading terms in equation (1) are singular like r ; all subsequent 
terms are nonsingular. The coefficients s^ and are related to the opening 

and sliding mode stress-intensity factors Kj and by the following formulas: 
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The strains corresponding to the stresses in equation (1) are obtained through 
Hooke's law. The strain-displacement equations in plane polar coordinates can 
then be integrated for the radial and tangential displacement components u^ (r,6) 

and Uq (r,0), respectively. See, for example, reference 7 and the subsequent 

correction. The resulting displacement components are 
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in which F, H, and K are rigid-body displacement parameters and G is the shear 
modulus. The dimensionless elastic constant a is given by 

(for plane strain) 

(for plane stress) 

where v is Poisson's ratio. 

Most cracked finite elements developed to date (c.f. refs. 8-11) incorpo- 
rate only the leading synmetric terms in equations (1) and (4). Creager, at the 
Lockheed-California Company in 1970, attempted to include subsequent terms in 
the Williams' series out was unsuccessful due to inadequate element geometry. 
This was successfully accomplished by Wilson (ref. 12) with a symmetric ele- 
ment, which makes use of the first four terms. Wilson's element, however, has 
the disadvantage of being semicircular and, hence, is awkward to use in con- 
junction with conventional elements which almost always have straight bound- 
aries. Moreover, the Wilson element (as well as some others previously 
referenced) has fewer degrees of freedom than are needed for independence of 
the nodal displacements. At best this requires that the stiffness matrix of 
the cracked element receive special attention in forming the stiffness matrix 
of the assembly. 

At Lockheed-Georgia, the decision was made at the outset to develop a 
cracked finite element that is a high-order element in that it 
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(i) incorporates many of the terms in the Williams' series; 

(ii) has a perfect balance between actual degrees of freedom and number 
of nodal displacements; and 

(iii) has a convenient shape for interfacing with conventional elements. 

The first feature permits accurate estimates of stress-intensity factors 
with relatively coarse finite-element grids, while the second and third fea- 
tures allow the numerical analyst to add the cracked element to an assembly 
in the same way that he adds a conventional element. 

Because many fracture mechanics problems are symmetric about the plane of 
the crack, two elements were developed at Lockheed-Georgia, One takes only 
the symmetric terms in the Williams' series and, hence, is applicable only to 
symmetric problems (Kji * 0) ; the other makes use of both symmetric and anti- 
symmetric terms and is applicable to unsymmetric or mixed mode (Kj and Kxi) 
problems. 


Plane-Deformation Symmetric Element 

Figure 2 shows the eight-node-syrametric element. The elemental coordinate 
system has its origin at the crack tip. The element is rectangular with a 
three-to-one aspect ratio. Placement of the nodes relative to the rectangle is 
pre-determined with a node at each corner plus nodes at the one-third points 
of each of the long sides. The choice of the symmetric element's shape - three 
equal squares - was considered convenient since the use of regular mesh spac- 
ings is common in finite-element models. This geometry has also proven to be 
effective when used with constant- strain triangles. The lower side (nodes 6, 

7, 8, and 1) is coincident with the crack direction and presumed axis of sym- 
metry. Nodes 6 and 7 are on the free crack face. Nodes 8 and 1 are on the 
prolongation of the crack. They are constrained rigidly as to vertical dis- 
placement and are free of shear forces - conditions consistent with symmetry. 

The element has sixteen displacement degrees of freedom - two per node 
corresponding to the in-plane displacement components. Thus, in keeping with 
feature (ii) mentioned previously, it incorporates the first thirteen s>mmet- 
ric terms of the Williams' series plus the three displacement degrees of free- 
dom associated with rigid-body displacement in the plane. In the following, 
the thirteen Williams' coefficients and the three rigid-body parameters are 
referred to as the sixteen generalized coordinates of the element. The 
stresses and displacements corresponding to these sixteen generalized coordin- 
ates are evaluated on the boundary of the element. Products of .'tress and dis- 
placement contributing to boundary work are formed and integrated. The result 
is a homogeneous quadratic form in the generalized coordinates, and the co- 
efficient of each term is an element of the cracked element's stiffness matrix 
with resp' t to the generalized coordinates. 

Once the stiffness matrix relative to generalized coordinates is deter- 
mined, the stiffness matrix relative to nodal displacements is formed using 
equation (4). 



To date this element has not failed to substantiate any reliable stress- 
intensity factor to within 2 % difference. And more often than not, this was 
accomplished with a relatively coarse finite-element grid. Results obtained 
with the symmetric element will be discussed later in the section "Applications." 


Plane-Deformation Unsymmetric Element 

The ten-node unsymmetric element is shown in figure 3. The element is 
square with equally spaced nodes around its boundary. As before, the shape 
(4 equal squares) and relative location of the nodes are fixed, and were 
chosen to provide modeling convenience. Again, the actual size of the element 
and its elastic constants dictated by the particular application are input 
parameters. The generalized coordinates correspond to the first nine symmet- 
ric terms and first eight antisymmetric terms of the Williams' series plus the 
three rigid-body displacement parameters. The stiffness matrix was again 
generated by integration around the boundary. 

Results obtained with this element are presented in the following section. 
Although sufficient, the accuracy obtained with this unsymmetric element is 
not as impressive as that obtained with the symmetric element previously dis- 
cussed. This is understandable in light of the fact that the unsymmetric ele- 
ment has fewer degrees of freedom that it can bring to bear on each mode. Of 
course, it can be used in a much wider class of crack problems and is more 
practical for industrial applications. In the following section some examples 
are given which show the accuracy and applicabilities of both the symmetric 
and unsymmetric elements. 


APPLICATIONS 


To illustrate the capabilities of the two elements, four examples of their 
usage follow. These examples were chosen to demonstrate first, the accuracy 
and economy of the elements and second, the capacities of the elements to per- 
form analyses for structural configurations of practical importance. The four 
cases are drawn from a wide range of work involving use of the cracked ele- 
ments and represent typical rather than most favorable results. 

The cracked elements are implemented in a standard single-precision 
finite-element displacement -method program which employs a banded Cholesky 
decomposition solution procedure. This program, which operates on a Univac 
1106 computer, was used for all four examples. 


Case 1: Symmetric- Element Test Case 

This example, shown on figure 4, was one of the first analyzed with the 
symmetric element. However, it exhibits the degree of accuracy which has been 
consistently achieved in numerous subsequent problems. The finite-element 
model has 31 nodes, 35 constant -strain triangles, and 1 eight-node cracked 
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element. The three configurations of single-edge crack, double-edge crack 
and center crack were all individually analyzed with this one model for an 
a/w ratio of 1/3. The model grid, which is quite coarse, results in the 
single-edge crack model having 57 displacement degrees of freedom (DOF) while 
the double-edge and center models have 51 DOF each. Comparisons of the 
stress-intensity factors computed using these conf igurations with ASTM values 
are shown on figure 4. The accuracy of the finite-element predictions are 
impressive (< 1.5% error) for all three cases. Computer time to perform each 
analysis was 3 to 4 seconds. 

Subsequent work with this finite-element model and others like it showed 
comparable results, with refinements in the grid bringing steady convergence 
toward ASTM values. 


Case 2: Symmetric Cracked Hole 

The symmetric problem depicted in figure 5 was analyzed in order to assess 
the accuracy of the formula 


v - Q 

I ~ H Bowie 


j3 T , , a ,/na 
Isida 


(5) 


which is routinely used to estimate stress-intensity factors for this geometry. 
In equation (5) B fiowie and 6 Isida are the correction factors, respectively, 

associated with Bowie’s (ref. 13) analysis for the presence of the hole and 
Isida' s (ref. 14) analysis for the finite-width effect. The finite-element 
grid as shown in figure 5 was established (i) to permit the location of the 
cracked element to be easily changed to simulate growth of the symmetric 
crack, and (ii) to permit the width of the plate to be readily changed by re- 
moving columns of constant-strain triangular elements from the right edge. The 
stress-intensity factors computed in this parametric study are listed in 
table 1. An important conclusion to be drawn from these results is that the 
value predicted by equation (5) seems adequate for short cracks, but becomes 
quite inaccurate and nonconserv&tive as the crack grows long. 

Such parametric analyses are of practical value, however, only if they 
can be quickly and economically accomplished. The use of the symmetric crack- 
ed element made it possible for a relatively coarse model (figure 5) to give 
accurate answers. The analysis of this coarse model required 28 000 words of 
computer storage and approximately 20 seconds of computer time. The total 
study, which includes model conception, data input, and 21 separate 
analyses, required less than 1 man-day of engineering effort and 7 minutes of 
computer time. Subsequent studies of other parameters, e.g, fastener loads 
or fastener interference, could now oe accomplished with even less effort 
since the model can be saved and reused. 

In general, the computer times and storage cited above are economical, 
and, in addition, are well within the operational limits of most remote access 
or time-sharing computer facilities. This economy and ready availability in a 
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remote access mode (Univac DEMAND) have resulted in widespread application of 
cracked elements at the Lockheed-Georgia Company in studies such as the one 
described above. 


Case 3: Pin-Loaded Lug 

The structural member shown on figure 6 is a lug and occurs frequently in 
aerospace structures The lug's lack of geometric symmetry greatly hampers 
the use of approximate methods to adequately estimate its stress-intensity 
factors Kj and Kxi when a crack appears. Such cracked geometry is easily 
handled, however, with the unsymmetric cracked element. The lug and the pin 
loading it are modeled with constant -strain triangles and a cracked element as 
shown on figure 7. The unsymmetric cracked element is shown in its initial 
position. As in the previous example, the model grid was constructed to per- 
mit the position of the cracked element to be easily changed to simulate the 
growth of the crack. Results of the analyses are given in Table 2. The small 
value of Kxi relative to Ki for the initial position indicates that the be- 
havior of the crack is primarily mode I or opening mode from the beginning. 
Figure 8 depicts the mode I behavior as the crack progresses from the hole to 
the outer edge. 

This analysis, which is considered to be in error by less than 3 percent, 
was accomplished at the cost of 1 man-day and 4 minutes of computer time in 
35 000 words of computer storage. 


Case 4: 45-Degree Slanted Crack Test Specimen 

This last example demonstrates the use of cracked elements to calibrate 
test configurations for accurate reduction of test data. The specimen shown 
by figure 9 contained an initial 45° center crack, and was subjected to a 
constant-amplitude tension-tension load. The maximum tension to minimum 
tension ratio was 0.1. The specimen was modeled with constant- strain trian- 
gles and an unsymmetric cracked element as shown by figure 10. The path taken 
by the crack during the cyclic test was simply traced on the finite-element 
model. The center section of the model was reconstructed seven times to ac- 
comodate the cracked elements for dimensionless projected crack lengths X = 
a/w of 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, and 0.8. The steel loading pins for the 
test were modeled by spring elements spread over the approximate bearing sur- 
faces as shown. 

The results from ths analyses of this model are shown on figure 11. The 
curves on figure 11 are considered accurate within 2 percent. These results 
show an interesting and potentially significant feature not found in similar 
work reported by Iida and Kobayashi (ref. 15). For 0.5 £ X s 0.75 the mode 1 
stress-intensity factor Ki is greater than Ki for a straight crack of the same 
projected length. This increase in Ki (the "hump" in figure 11) occurs when 
Kn goes to zero and could easily account for the increased crack-growth rates 
observed for this type of specimen. However, the significant point to this 
analysis is that an accurate knowledge of Kj and Kji for test articles permits 
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a concise reduction of crack -growth- rate data for later design usage. 

The analysis of this case required approximately 80 seconds of computer 
time for each configuration. The total for all seven configurations, including 
data input and modifications, was 9.4 minutes. 


CONCLUSIONS 


The incorporation into NASTRAN and other finite element codes of cracked 
elements appears to be both a timely and practical effort, for there is a grow- 
ing concern among those involved in aerospace design with being able to per- 
form reliable fracture analyses of damaged or flawed structures. This concern 
arises partly from anticipation of difficulties of meeting certain design 
"fracture criteria" imposed by procuring agencies for new aerospace vehicles 
and partly from recent experiences with existing airplanes. For whatever 
reasons, the fact is that fracture mechanics and fracture analyses have become 
significant and necessary steps in the design and modifications of aerospace 
structures. Thus, there is a pressing need for accurate production analysis 
tools to enable designers to apply new fracture criteria. Cracked finite 
elements, such as those discussed in this paper, are such tools. Their ease 
of application makes it possible for analysts and designers not extensively 
trained in fracture mechanics but familiar with finite element methods to com- 
pute accurate stress-intensity factors. 
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Table 1. A Study of Finite-Width Panel with Two Symmetric Cracks Coming From a Hole 
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Table 2. Results from Analysis of the Lug Model 
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Figure 4.- Results for single-edge, double-edge, and center cracked tension panels 

with a/w = 1/3. 
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Figure 5.- Symmetric cracks growing from a hole in a finite -width panel. 
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Figure 10.- Finite -element models of tension specimen where initial 45 crack grew 

experimentally as modeled. 









